/****************************************************************************
 * Copyright (C) 2009-2010 SciTouch LLC
 * 
 * This file is part of Indigo toolkit.
 * 
 * This file may be distributed and/or modified under the terms of the
 * GNU General Public License version 3 as published by the Free Software
 * Foundation and appearing in the file LICENSE.GPL included in the
 * packaging of this file.
 * 
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 ***************************************************************************/

#include <stdarg.h>
#include <stdlib.h>
#include <stdio.h>

#include "base_c/defs.h"
#include "base_cpp/tlscont.h"

#include "molecule/molecule.h"
#include "math/algebra.h"
#include "graph/spanning_tree.h"
#include "molecule/molecule_stereocenters.h"
#include "molecule/molecule_3d_constraints.h"
#include "graph/filter.h"

Molecule::Molecule ()
{
   _stereocenters.create(*this);
   _3d_constraints = 0;
   have_xyz = false;
}

Molecule::~Molecule ()
{
   delete _3d_constraints;
}

void Molecule::clear ()
{
   Graph::clear();

   _atoms.clear();
   _bonds.clear();
   _connectivity.clear();
   _bonds_in_rings.clear();
   _atoms_in_aromatic_rings.clear();
   _stereocenters->clear();
   delete _3d_constraints;
   _3d_constraints = 0;
   _rgroups.reset(0);
   _rgroup_fragment.reset(0);
   fixed_atoms.clear();
   cis_trans.clear();
   if (_q_atoms.get() != 0)
      _q_atoms->clear();
   if (_q_bonds.get() != 0)
      _q_bonds->clear();
   have_xyz = false;
   _pseudo_atom_values.clear();
   _pseudo_atoms.clear();
   name.clear();
}

void Molecule::registerAtom (int idx, const Atom &atom)
{
   _atoms.expand(idx + 1);

   _atoms[idx] = atom;

   _atoms_in_aromatic_rings.expand(idx + 1);
   _atoms_in_aromatic_rings[idx] = 0;
   if (_q_atoms.get() != 0)
   {
      _q_atoms->expand(idx + 1);
      _q_atoms->at(idx).clear();
   }
}

void Molecule::registerBond (int idx, const Bond &bond)
{
   _bonds.expand(idx + 1);
   _bonds[idx] = bond;

   _bonds_in_rings.expand(idx + 1);
   _bonds_in_rings[idx] = 0;

   if (bond.type == BOND_AROMATIC)
      setBondAromatic(idx);

   cis_trans.registerBond(idx);
   
   if (_q_bonds.get() != 0)
   {
      _q_bonds->expand(idx + 1);
      _q_bonds->at(idx).clear();
   }
}

int Molecule::addAtom (int label)
{
   Atom new_atom;

   new_atom.label = label;

   return addAtomCopy(new_atom);
}

int Molecule::addAtomCopy (const Atom &atom)
{
   int res = addVertex();

   registerAtom(res, atom);

   return res;
}

int Molecule::addBond (int beg, int end, int type)
{
   return addBondCopy(beg, end, Bond(type));
}

int Molecule::addBondCopy (int beg, int end, const Bond &bond)
{
   int res = addEdge(beg, end);

   registerBond(res, bond);

   return res;
}

void Molecule::flipBond (int atom_parent, int atom_from, int atom_to)
{
   int src_bond_idx = findEdgeIndex(atom_parent, atom_from);
   Bond src_bond = getBond(src_bond_idx);
   int added_idx = addBondCopy(atom_parent, atom_to, src_bond);
   if (haveQueryBonds())
      getQueryBond2(added_idx) = getQueryBond(src_bond_idx);

   getStereocenters().flipBond(atom_parent, atom_from, atom_to);
   cis_trans.flipBond(*this, atom_parent, atom_from, atom_to);
   removeEdge(src_bond_idx);
}

void Molecule::_mergeWithSubmolecule (const Molecule &molecule, const Array<int> &vertices,
                                      const Array<int> *edges, Array<int> *mapping_out)
{
   QS_DEF(Array<int>, tmp_mapping);

   if (mapping_out == 0)
      mapping_out = &tmp_mapping;

   int i;

   if (vertexCount() == 0)
      have_xyz = molecule.have_xyz;
   else
      have_xyz = have_xyz && molecule.have_xyz;

   // vertices and edges
   _mergeWithSubgraph(molecule, vertices, edges, mapping_out);

   if (molecule.haveQueryAtoms())
      createQueryAtoms();
   if (molecule.haveQueryBonds())
      createQueryBonds();

   // atoms
   for (i = 0; i < vertices.size(); i++)
      registerAtom(mapping_out->at(vertices[i]), molecule.getAtom(vertices[i]));

   // pseudo-atoms
   for (i = 0; i < vertices.size(); i++)
   {
      if (molecule.isPseudoAtom(vertices[i]))
         setPseudoAtom(i, molecule.getPseudoAtom(vertices[i]));
   }

   // query atoms
   if (molecule.haveQueryAtoms())
      for (i = 0; i < vertices.size(); i++)
         _q_atoms->at(mapping_out->at(vertices[i])) = molecule.getQueryAtom(vertices[i]);

   // bonds and query bonds and bonds in rings
   if (edges != 0)
      for (i = 0; i < edges->size(); i++)
      {
         const Edge &edge = molecule.getEdge(edges->at(i));
         int beg = mapping_out->at(edge.beg);
         int end = mapping_out->at(edge.end);

         if (beg == -1 || end == -1)
            // must be thrown before in mergeWithSubgraph()
            throw Error("_mergeWithSubmolecule: internal");

         int idx = findEdgeIndex(beg, end);

         registerBond(idx, molecule.getBond(edges->at(i)));

         if (molecule.haveQueryBonds())
            _q_bonds->at(idx) = molecule.getQueryBond(edges->at(i));

         if (molecule.bondInRing(edges->at(i)))
            _bonds_in_rings[idx] = 1;
      }
   else for (i = molecule.edgeBegin(); i < molecule.edgeEnd(); i = molecule.edgeNext(i))
   {
      const Edge &edge = molecule.getEdge(i);
      int beg = mapping_out->at(edge.beg);
      int end = mapping_out->at(edge.end);

      if (beg == -1 || end == -1)
         continue;

      int idx = findEdgeIndex(beg, end);

      registerBond(idx, molecule.getBond(i));

      if (molecule.haveQueryBonds())
         _q_bonds->at(idx) = molecule.getQueryBond(i);

      if (molecule.bondInRing(i))
         _bonds_in_rings[idx] = 1;
   }

   // aromatic rings
   _atoms_in_aromatic_rings.clear_resize(vertexEnd());
   _atoms_in_aromatic_rings.zerofill();

   for (i = 0; i < vertices.size(); i++)
      if (molecule.atomInAromaticRing(vertices[i]))
         _atoms_in_aromatic_rings[mapping_out->at(vertices[i])] = 1;

   cis_trans.buildOnSubmolecule(molecule, *this, mapping_out->ptr());

   _stereocenters->buildOnSubmolecule(molecule._stereocenters.ref(), mapping_out->ptr());

   // 3d constraints
   if (molecule._3d_constraints != 0)
   {
      Molecule3dConstraints *constraints = create3dConstraints();

      constraints->buildOnSubmolecule(*molecule._3d_constraints, mapping_out->ptr());

      if (!constraints->haveConstraints())
      {
         delete constraints;
         _3d_constraints = 0;
      }
   }

   // fixed atoms
   for (i = 0; i < molecule.fixed_atoms.size(); i++)
   {
      int idx = mapping_out->at(molecule.fixed_atoms[i]);

      if (idx >= 0)
         fixed_atoms.push(idx);
   }

   // RGroups and RGroup fragments
   if (molecule.getRGroupFragment() != 0)
   {
      MoleculeRGroupFragment &other_fragment = *molecule.getRGroupFragment();
      MoleculeRGroupFragment *fragment = getRGroupFragment();
      
      if (fragment == 0)
         fragment = createRGroupFragment();

      for (i = 0; i < other_fragment.attachmentPointCount(); i++)
      {
         int att_idx;
         int j;
         
         for (j = 0; (att_idx = other_fragment.getAttachmentPoint(i, j)) != -1; j++)
            fragment->addAttachmentPoint(i, mapping_out->at(att_idx));
      }
   }

   if (molecule.getRGroups() != 0)
   {
      MoleculeRGroups *rgroups = getRGroups();

      if (rgroups == 0)
         rgroups = createRGroups();

      rgroups->copyRGroupsFromMolecule(*molecule.getRGroups());
      rgroups->copySitesFromMolecule(*molecule.getRGroups(), mapping_out->ptr());
   }
}


void Molecule::makeSubmolecule (const Molecule &molecule,
                                const Array<int> &vertices, Array<int> *mapping_out)
{
   clear();
   _mergeWithSubmolecule(molecule, vertices, 0, mapping_out);
}

void Molecule::makeSubmolecule (const Molecule &other, const Filter &filter,
                                Array<int> *mapping_out, Array<int> *inv_mapping)
{
   QS_DEF(Array<int>, vertices);

   if (mapping_out == 0)
      mapping_out = &vertices;

   filter.collectGraphVertices(other, *mapping_out);

   makeSubmolecule(other, *mapping_out, inv_mapping);
}


void Molecule::makeEdgeSubmolecule (const Molecule &molecule, const Array<int> &vertices, const Array<int> &edges, Array<int> *v_mapping)
{
   clear();
   _mergeWithSubmolecule(molecule, vertices, &edges, v_mapping);
}

void Molecule::cloneMolecule (const Molecule &other, Array<int> *mapping, Array<int> *inv_mapping)
{
   QS_DEF(Array<int>, tmp_mapping);

   if (mapping == 0)
      mapping = &tmp_mapping;

   mapping->clear();

   for (int i = other.vertexBegin(); i < other.vertexEnd(); i = other.vertexNext(i))
      mapping->push(i);

   makeSubmolecule(other, *mapping, inv_mapping); 
}  

void Molecule::setQueryBondAromatic (int idx)
{
   setBondAromatic(idx);
   _q_bonds->at(idx).type = 0;
   _q_bonds->at(idx).can_be_aromatic = false;
}

void Molecule::setQueryBondFuzzyAromatic (int idx)
{
   _q_bonds->at(idx).can_be_aromatic = true;
}

const Atom & Molecule::getAtom (int idx) const
{
   return _atoms[idx];
}

Atom & Molecule::getAtom2 (int idx)
{
   return _atoms[idx];
}

const Bond & Molecule::getBond (int idx) const
{
   return _bonds[idx];
}

Bond & Molecule::getBond2 (int idx)
{
   return _bonds[idx];
}

void Molecule::setBondType (int idx, int type)
{
   _bonds[idx].type = type;
}

void Molecule::setAtom (int idx, const Atom &atom)
{
   _atoms[idx] = atom;
}

void Molecule::setAtomCharge (int idx, int charge)
{
   _atoms[idx].charge = charge;
}

const Vec3f & Molecule::getAtomPos (int idx) const
{
   return _atoms[idx].pos;
}

void Molecule::setAtomPos (int idx, const Vec3f& new_pos)
{
   _atoms[idx].pos = new_pos;
}

void Molecule::setPseudoAtom (int idx, const char *text)
{
   _atoms[idx].label = 0;
   if (_pseudo_atoms.at2(idx) != 0)
      _pseudo_atoms.at(idx) = _pseudo_atom_values.add(text);
   else
      _pseudo_atoms.insert(idx, _pseudo_atom_values.add(text));
}

bool Molecule::isPseudoAtom (int idx) const
{
   return _pseudo_atoms.at2(idx) != 0;
}

const char * Molecule::getPseudoAtom (int idx) const
{
   return _pseudo_atom_values.at(_pseudo_atoms.at(idx));
}


void Molecule::mergeWithMolecule (const Molecule &other, Array<int> *mapping)
{
   QS_DEF(Array<int>, vertices);
   int i;

   vertices.clear();

   for (i = other.vertexBegin(); i != other.vertexEnd(); i = other.vertexNext(i))
      vertices.push(i);

   _mergeWithSubmolecule(other, vertices, 0, mapping);
}

int Molecule::getVacantPiOrbitals (int atom_idx, int conn, int *lonepairs_out) const
{
   const Atom &atom = getAtom(atom_idx);
   int group = Element::group(atom.label);
   int charge = atom.charge;

   int orbitals;
   switch (group)
   {
      case 1: orbitals = 1; break;
      case 2: orbitals = 2; break;
      default: orbitals = 4;
   }

   int free_electrons = group - conn - charge - atom.radical;
   if (free_electrons < 0)
      return -1;

   int lonepair = free_electrons / 2;
   int implicit_radical = free_electrons % 2;

   int vacant = orbitals - conn - lonepair - atom.radical - implicit_radical;
   if (vacant < 0)
      return -1;

   if (lonepairs_out != 0)
      *lonepairs_out = lonepair;
   return vacant;
}

int Molecule::getVacantPiOrbitals (int atom_idx, int *lonepairs_out) const
{
   return getVacantPiOrbitals(atom_idx, _connectivity[atom_idx], lonepairs_out);
}

int Molecule::totalCharge () const
{
   int i, total_charge = 0;

   for (i = vertexBegin(); i < vertexEnd(); i = vertexNext(i))
      total_charge += getAtom(i).charge;

   return total_charge;
}

int Molecule::totalHydrogensCount () const
{
   int i, total_h = 0;

   for (i = vertexBegin(); i < vertexEnd(); i = vertexNext(i))
   {
      const Atom &atom = getAtom(i);

      if (atom.label == ELEM_H)
         total_h++;

      total_h += atom.implicit_h;
   }

   return total_h;
}

bool Molecule::matchAtomsEq (const Graph &g1, const Graph &g2,
                             int idx1, int idx2, void *userdata)
{
   return matchAtomsCmp(g1, g2, idx1, idx2, userdata) == 0;
}

int Molecule::matchAtomsCmp (const Graph &g1, const Graph &g2,
                             int idx1, int idx2, void *userdata)
{
   const Molecule &m1 = (const Molecule &)g1;
   const Molecule &m2 = (const Molecule &)g2;

   const Atom &atom1 = m1.getAtom(idx1);
   const Atom &atom2 = m2.getAtom(idx2);

   if (m1.isPseudoAtom(idx1) && !m2.isPseudoAtom(idx2))
      return 1;

   if (!m1.isPseudoAtom(idx1) && m2.isPseudoAtom(idx2))
      return -1;

   if (m1.isPseudoAtom(idx1) && m2.isPseudoAtom(idx2))
   {
      int res = strcmp(m1.getPseudoAtom(idx1), m2.getPseudoAtom(idx2));

      if (res != 0)
         return res;
   }
   else
   {
      if (atom1.label > atom2.label)
         return 1;
      if (atom1.label < atom2.label)
         return -1;
   }

   if (atom1.isotope > atom2.isotope)
      return 1;
   if (atom1.isotope < atom2.isotope)
      return -1;

   if (atom1.charge > atom2.charge)
      return 1;
   if (atom1.charge < atom2.charge)
      return -1;

   if (atom1.valence > atom2.valence)
      return 1;
   if (atom1.valence < atom2.valence)
      return -1;

   if (atom1.radical > atom2.radical)
      return 1;
   if (atom1.radical < atom2.radical)
      return -1;

   return 0;
}

bool Molecule::matchBondsEq (const Graph &g1, const Graph &g2,
                              int idx1, int idx2, void *userdata)
{
   const Molecule &m1 = (const Molecule &)g1;
   const Molecule &m2 = (const Molecule &)g2;

   const Bond &bond1 = m1.getBond(idx1);
   const Bond &bond2 = m2.getBond(idx2);

   if (bond1.type != bond2.type)
      return false;

   return true;
}

void Molecule::setBondAromatic (int idx)
{
   const Edge &edge = _edges[idx];
   Bond &bond = _bonds[idx];

   bond.type = BOND_AROMATIC;

   _atoms_in_aromatic_rings[edge.beg] = 1;
   _atoms_in_aromatic_rings[edge.end] = 1;
}

bool Molecule::atomInAromaticRing (int idx) const
{
   return _atoms_in_aromatic_rings[idx] != 0;
}

bool Molecule::bondInRing (int idx) const
{
   return _bonds_in_rings[idx] != 0;
}

void Molecule::setBondInRing (int idx, bool in_ring)
{
   _bonds_in_rings[idx] = in_ring ? 1 : 0;
}

void Molecule::makeSubmoleculeFromQuery (const Molecule &query, Array<int> *mapping_out, Array<int> *inv_mapping_out)
{
   QS_DEF(Array<int>, mapping_tmp);
   QS_DEF(Array<int>, inv_mapping_tmp);
   int i;

   if (mapping_out == 0)
      mapping_out = &mapping_tmp;
   if (inv_mapping_out == 0)
      inv_mapping_out = &inv_mapping_tmp;

   mapping_out->clear();

   for (i = query.vertexBegin(); i < query.vertexEnd(); i = query.vertexNext(i))
      if (query.getQueryAtom(i).type == 0)
         mapping_out->push(i);

   makeSubmolecule(query, *mapping_out, inv_mapping_out);

   // remove the query bonds from the subgraph
   for (i = query.edgeBegin(); i < query.edgeEnd(); i = query.edgeNext(i))
   {
      const Edge &edge = query.getEdge(i);

      if (query.getQueryBond(i).type == 0)
         continue;

      int idx1 = inv_mapping_out->at(edge.beg);
      int idx2 = inv_mapping_out->at(edge.end);

      if (idx1 == -1 || idx2 == -1)
         continue; // was removed with vertex

      removeEdge(findEdgeIndex(idx1, idx2));
   }
}

void Molecule::findBondsInRings ()
{
   SpanningTree spt(*this);

   _bonds_in_rings.zerofill();
   spt.markAllEdgesInCycles(_bonds_in_rings.ptr(), 1);
}

bool Molecule::haveZCoord (const Molecule &mol)
{
   for (int i = mol.vertexBegin(); i != mol.vertexEnd(); i = mol.vertexNext(i))
      if (fabs(mol.getAtom(i).pos.z) > 0.001)
         return true;

   return false;
}

RGroup::RGroup () : if_then(-1), rest_h(0)
{
}

RGroup::~RGroup ()
{
}

void RGroup::clear()
{
   if_then = -1;
   rest_h = 0;
   occurence.clear();
   fragments.clear();
}

MoleculeStereocenters & Molecule::getStereocenters () const
{
   return _stereocenters.ref();
}

Molecule3dConstraints * Molecule::create3dConstraints ()
{
   if (_3d_constraints != 0)
      throw Error("attempt to create 3d constraints twice");

   _3d_constraints = new Molecule3dConstraints(*this);
   return _3d_constraints;
}

Molecule3dConstraints * Molecule::get3dConstraints () const
{
   return _3d_constraints;
}

MoleculeRGroups * Molecule::createRGroups()
{
   if (_rgroups.get() != 0)
      throw Error("attempt to create rgroups twice");

   _rgroups.reset(new MoleculeRGroups(*this));
   return _rgroups.get();
}

MoleculeRGroups * Molecule::getRGroups() const
{
   return _rgroups.get();
}

MoleculeRGroupFragment * Molecule::createRGroupFragment()
{
   if (_rgroup_fragment.get() != 0)
      throw Error("attempt to create rgroup fragment twice");

   _rgroup_fragment.reset(new MoleculeRGroupFragment());
   return _rgroup_fragment.get();
}

MoleculeRGroupFragment * Molecule::getRGroupFragment() const
{
   return _rgroup_fragment.get();
}

void Molecule::clearRGroupFragment()
{
   _rgroup_fragment.reset(0);
}

void Molecule::removeAtoms (const Array<int> &indices)
{
   QS_DEF(Array<int>, mapping);
   int i;

   _stereocenters->removeAtoms(indices);

   mapping.clear_resize(vertexEnd());

   for (i = vertexBegin(); i != vertexEnd(); i = vertexNext(i))
      mapping[i] = i;

   for (i = 0; i < indices.size(); i++)
   {
      mapping[indices[i]] = -1;
      removeVertex(indices[i]);
   }

   if (_3d_constraints != 0)
   {
      AutoPtr<Molecule3dConstraints> old(_3d_constraints);

      _3d_constraints = 0; // to keep it clear even if the constructor call throws an error
      _3d_constraints = new Molecule3dConstraints(*this);
      _3d_constraints->buildOnSubmolecule(old.ref(), mapping.ptr());
   }

   cis_trans.buildOnSubmolecule(*this, *this, mapping.ptr());

   if (_rgroups.get() != 0)
      _rgroups->removeSites(indices);

   if (_rgroup_fragment.get() != 0)
   {
      for (i = 0; i < indices.size(); i++)
         _rgroup_fragment->removeAttachmentPoint(indices[i]);

      bool empty = true;

      for (i = 0; i < _rgroup_fragment->attachmentPointCount(); i++)
         if (_rgroup_fragment->getAttachmentPoint(i, 0) != -1)
         {
            empty = false;
            break;
         }

      if (empty)
         clearRGroupFragment();
   }
}

void Molecule::removeAtom (int idx)
{
   QS_DEF(Array<int>, vertices);

   vertices.clear();
   vertices.push(idx);
   removeAtoms(vertices);
}

void Molecule::removeAtoms (const Filter &filter)
{
   QS_DEF(Array<int>, vertices);

   filter.collectGraphVertices(*this, vertices);
   removeAtoms(vertices);
}

void Molecule::removeBonds (const Array<int> &indices)
{
   _stereocenters->removeBonds(indices);

   for (int i = 0; i < indices.size(); i++)
      removeEdge(indices[i]);
}

void Molecule::removeBond (int idx)
{
   QS_DEF(Array<int>, edges);

   edges.clear();
   edges.push(idx);
   removeBonds(edges);
}

void Molecule::unfoldHydrogens (Array<int> *markers_out)
{
   int i, j, v_end = vertexEnd();

   if (markers_out != 0)
   {
      markers_out->clear_resize(vertexEnd());
      markers_out->zerofill();
   }

   for (i = vertexBegin(); i < v_end; i = vertexNext(i))
   {
      for (j = 0; j < _atoms[i].implicit_h; j++)
      {
         int new_h_idx = addAtom(ELEM_H);

         addBond(i, new_h_idx, BOND_SINGLE);
         if (markers_out != 0)
         {
            markers_out->expand(new_h_idx + 1);
            markers_out->at(new_h_idx) = 1;
         }

         getStereocenters().registerUnfoldedHydrogen(i, new_h_idx);
      }

      _atoms[i].implicit_h = 0;
   }
}

void Molecule::createQueryAtoms ()
{
   if (_q_atoms.get() != 0)
      return;
      
   _q_atoms.create();
   _q_atoms->clear_resize(vertexEnd());
   _q_atoms->zerofill();
}

void Molecule::createQueryBonds ()
{
   if (_q_bonds.get() != 0)
      return;
   
  _q_bonds.create();
  _q_bonds->clear_resize(edgeEnd());
  _q_bonds->zerofill();
}

const QueryAtom & Molecule::getQueryAtom (int idx) const
{
   if (_q_atoms.get() == 0)
      throw Error("attempt to get query atom of non-query molecule");
              
   return _q_atoms->at(idx);
}

const QueryBond & Molecule::getQueryBond (int idx) const
{
   if (_q_bonds.get() == 0)
      throw Error("attempt to get query bond of non-query molecule");
              
   return _q_bonds->at(idx);
}

QueryAtom & Molecule::getQueryAtom2 (int idx)
{
   if (_q_atoms.get() == 0)
      throw Error("attempt to get query atom of non-query molecule");
              
   return _q_atoms->at(idx);
}

QueryBond & Molecule::getQueryBond2 (int idx)
{
   if (_q_bonds.get() == 0)
      throw Error("attempt to get query bond of non-query molecule");
              
   return _q_bonds->at(idx);
}

bool Molecule::haveQueryAtoms () const
{
   return _q_atoms.get() != 0;
}

bool Molecule::haveQueryBonds () const
{
   return _q_bonds.get() != 0;
}

void Molecule::calcConnectivity (bool to_throw)
{
   int i;
   
   _connectivity.clear_resize(vertexEnd());
   _connectivity.zerofill();
   
   for (i = edgeBegin(); i != edgeEnd(); i = edgeNext(i))
   {
      const Edge &edge = getEdge(i);
      const Bond &bond = getBond(i);
      int cardinality = 0;
      
      if (bond.type == BOND_SINGLE)
         cardinality = 1;
      else if (bond.type == BOND_DOUBLE)
         cardinality = 2;
      else if (bond.type == BOND_TRIPLE)
         cardinality = 3;
      else if (!to_throw)
         cardinality = -1;
      else
         throw Error("calcConnectivity(): bad bond type: %d", bond.type);
         
      if (cardinality >= 0)
      {
         if (_connectivity[edge.beg] >= 0)
            _connectivity[edge.beg] += cardinality;
         if (_connectivity[edge.end] >= 0)
            _connectivity[edge.end] += cardinality;
      }
      else
      {
         _connectivity[edge.beg] = -1;
         _connectivity[edge.end] = -1;
      }
   }
}

void Molecule::calcImplicitHydrogens (bool calc_maximum)
{
   int i;
   
   calcConnectivity(!calc_maximum);
   
   for (i = vertexBegin(); i != vertexEnd(); i = vertexNext(i))
   {
      Atom &atom = getAtom2(i);
      int conn = _connectivity[i];

      if (conn < 0)
      {
         atom.implicit_h = -1;
         continue;
      }

      if (atom.explicit_valence)
      {
         atom.implicit_h = atom.valence - Element::calcValenceMinusHyd(atom.label, atom.charge, atom.radical, conn);
         
         if (atom.implicit_h < 0)
         {
            if (calc_maximum)
            {
               atom.implicit_h = -1;
               continue;
            }
               
            throw Error("valence %d specified on %s, charge %d, radical %d, but %d bonds are drawn",
                    atom.valence, Element::toString(atom.label), atom.charge, atom.radical, conn);
         }
      }
      else
         Element::calcValence(atom.label, atom.charge, atom.radical, conn, atom.valence, atom.implicit_h, !calc_maximum);
     
      _connectivity[i] += atom.implicit_h;
   }
} 

void Molecule::calcExplicitValencesByImplicitHydrogens ()
{
   int i;
   calcConnectivity(false);

   for (i = vertexBegin(); i != vertexEnd(); i = vertexNext(i))
   {
      const Atom &atom = getAtom(i);
      int conn = _connectivity[i];

      if (atom.implicit_h < 0)
         continue;

      if (conn < 0)
         continue;

      int normal_val, normal_hyd;

      Element::calcValence(atom.label, atom.charge, atom.radical, conn, normal_val, normal_hyd, false);

      if (atom.implicit_h != normal_hyd && atom.radical == 0)
      {
         // try to put additional radicals

         // first try a singlet radical
         if (Element::calcValence(atom.label, atom.charge, RADICAL_SINGLET,
                                  conn, normal_val, normal_hyd, false) &&
             atom.implicit_h == normal_hyd)
            _atoms[i].radical = RADICAL_SINGLET;
         else
         {
            // try a douplet radical
            if (Element::calcValence(atom.label, atom.charge, RADICAL_DOUPLET,
                                     conn, normal_val, normal_hyd, false) &&
                    atom.implicit_h == normal_hyd)
               _atoms[i].radical = RADICAL_DOUPLET;
         }
      }

      if (atom.implicit_h != normal_hyd)
      {
         // radicals have not help, set explicit valence
         int val = Element::calcValenceMinusHyd(atom.label, atom.charge, atom.radical, conn) + atom.implicit_h;

         if (val >= 0)
         {
            _atoms[i].valence = val;
            _atoms[i].explicit_valence = true;
         }
         
      }

   }
}

bool Molecule::hasStereoCareBonds (const Molecule &query)
{
   int i;

   for (i = query.edgeBegin(); i != query.edgeEnd(); i = query.edgeNext(i))
      if (query.getQueryBond(i).stereo_care)
         return true;

   return false;
}

bool Molecule::hasExplicitValences (const Molecule &query)
{
   int i;

   for (i = query.vertexBegin(); i != query.vertexEnd(); i = query.vertexNext(i))
      if (query.getAtom(i).explicit_valence)
         return true;

   return false;
}

int Molecule::getStereocentersOffset ()
{
   char dummy[sizeof(Molecule)];

   return (int)((char *)(&((Molecule *)dummy)->_stereocenters) - dummy);
}

void Molecule::saveBondTypes (const Molecule &mol, Array<int> &types)
{
   int i;

   types.clear_resize(mol.edgeEnd());
   for (i = mol.edgeBegin(); i != mol.edgeEnd(); i = mol.edgeNext(i))
      types[i] = mol.getBond(i).type;
}

void Molecule::loadBondTypes (Molecule &mol, Array<int> &types)
{
   int i;

   for (i = mol.edgeBegin(); i != mol.edgeEnd(); i = mol.edgeNext(i))
      mol.getBond2(i).type = types[i];
}

bool Molecule::fragmentHasQueryBonds (const Molecule &mol, const Filter *filter)
{
   if (!mol.haveQueryBonds())
      return false;

   for (int i = mol.edgeBegin(); i != mol.edgeEnd(); i = mol.edgeNext(i))
   {
      if (filter != 0 && !filter->valid(i))
         continue;

      const QueryBond &qbond = mol.getQueryBond(i);

      if (qbond.type != 0 || qbond.can_be_aromatic)
         return true;
   }

   return false;
}

bool Molecule::fragmentHasQueryBonds (const Molecule &mol, const Array<int> &edges)
{
   if (!mol.haveQueryBonds())
      return false;

   for (int i = 0; i < edges.size(); i++)
   {
      const QueryBond &qbond = mol.getQueryBond(edges[i]);

      if (qbond.type != 0 || qbond.can_be_aromatic)
         return true;
   }

   return false;
}

bool Molecule::fragmentHasQueryAtoms (const Molecule &mol, const Filter *filter)
{
   if (!mol.haveQueryAtoms())
      return false;

   for (int i = mol.vertexBegin(); i != mol.vertexEnd(); i = mol.vertexNext(i))
   {
      if (filter != 0 && !filter->valid(i))
         continue;

      if (mol.getQueryAtom(i).type != 0)
         return true;
   }

   return false;
}

bool Molecule::fragmentHasQueryAtoms (const Molecule &mol, const Array<int> &vertices)
{
   if (!mol.haveQueryAtoms())
      return false;

   for (int i = 0; i != vertices.size(); i++)
      if (mol.getQueryAtom(vertices[i]).type != 0)
         return true;

   return false;
}
